%pylab inline
import glob
from matplotlib import image
from clawpack.visclaw.JSAnimation import IPython_display
from matplotlib import animation
def init():
im.set_data(image.imread(filenames[0]))
return im,
def animate(i):
image_i=image.imread(filenames[i])
im.set_data(image_i)
return im,
Populating the interactive namespace from numpy and matplotlib
In 1955, Lighthill, Whitham and Richards proposed the a model, now known as LWR model for traffic flow problem. This model is a macroscopic traffic flow model.
$$\rho_t+f(\rho)_x=0,$$where $$f(\rho)=\rho v(\rho).$$
From the empirical data collected from freeway, the velocity can be modeled as:
$$v(\rho)=\left\{\begin{array}{lll} 1,&\text{if }0\leq\rho<\rho_m,\\ \gamma(\frac{1}{\rho}-1),&\text{if }\rho_m\leq\rho\leq 1. \end{array}\right.$$Hence
$$f(\rho)=\left\{\begin{array}{lll} \rho,&\text{if }0\leq\rho<\rho_m,\\ \gamma(1-\rho),&\text{if }\rho_m\leq\rho\leq 1. \end{array}\right.$$For this problem we use $\rho_m=0.5$ and $\gamma=0.5$.
rhom = 0.5
gamma = 0.5
def flux_freeway_origin(rho):
n = rho.size
f = zeros(n)
for i in range(n):
if rho[i] >= 0. and rho[i] < rhom:
f[i] = rho[i]
elif rho[i] >= rhom and rho[i] <= 1.:
f[i] = gamma*(1.-rho[i])
else:
print "Invalid density value!"
return f
rho1 = linspace(0.,rhom-0.01,50)
rho2 = linspace(rhom,1.,51)
figure(figsize=(10,6))
plot(rho1,flux_freeway_origin(rho1),'b')
plot(rho2,flux_freeway_origin(rho2),'b')
plot([rhom,rhom],[rhom,gamma*(1.-rhom)],'--b')
ylim(0.,0.5)
xlabel('rho')
ylabel('$f$')
title('Original flux fucntion')
<matplotlib.text.Text at 0x106291ed0>
The discontinuity can result the invalidation of the technique we learned from the class, since the mehod involves the deriavtive of the flux function. One method to redeem this is to use mollifier smooth the flux function. We use the mollifier as follows:
$$\eta_\epsilon(s)=\frac{1}{\epsilon}\eta\left(\frac{s}{\epsilon}\right),$$where
$$\eta(x)=\left\{\begin{array}{lll} C\exp\left(\frac{1}{x^2-1}\right)&\text{if }-1< x < 1,\\ 0&\text{otherwise},\end{array}\right.$$Here $C$ is a constant, $C\approx2.2522836.$
The modified flux function is the convolution of the mollifier and original function:
$$f_\epsilon=\eta_\epsilon\ast f,$$so
$$f_\epsilon(\rho)=\rho+(\gamma(1-\rho)-\rho)\int_{\rho_m-\epsilon}^\rho\eta_\epsilon(s-\rho_m)ds,$$with $0<\epsilon\ll1$.
It has been proven that, if $\eta_\epsilon(x)\rightarrow\delta(x)$ as $\epsilon\rightarrow 0$, then $f_\epsilon\rightarrow f$.
In the project, we select $\epsilon=0.01$.
def eta(x,eps):
C = 2.25228362104358101049978125556
eta = where(abs(x)<eps,C*exp(eps**2/(x**2-eps**2))/eps,0.)
return eta
def flux_freeway_modified(rho,eps):
from scipy.integrate import quad
n = rho.size
f = zeros(n)
for i in range(n):
# N = 500*(int((rho[i]-rhom+eps)/eps) + 1) # mesh cell for integration
# dx = (rho[i]-rhom+eps)/N # mesh size for integration
# intl = 0.
# for j in range(N):
# intl = intl + dx*eta((j+0.5)*dx-eps,eps)
intl = quad(eta,-eps,rho[i]-rhom,args=eps)[0]
f[i] = rho[i] + (gamma*(1.-rho[i])-rho[i])*intl
return f
%%system
mkdir _plots_ini_freeway
['mkdir: _plots_ini_freeway: File exists']
rho = linspace(0.,1.,202)
i = 0
for epspic in linspace(0.5,0.01,50):
clf()
plot(rho,flux_freeway_modified(rho,epspic),'b')
ylim(0.,0.5)
xlabel('rho')
ylabel('$f_\epsilon$')
title('$f_\epsilon$ when $\epsilon$=%01f'%epspic)
fname = '_plots_ini_freeway/frame%04i'%i + 'fig1' + '.png'
savefig(fname)
i = i + 1
/Users/Peng/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: divide by zero encountered in double_scalars app.launch_new_instance() /Users/Peng/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: overflow encountered in exp app.launch_new_instance() /Users/Peng/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: overflow encountered in double_scalars app.launch_new_instance()
figno = 1
fname = '_plots_ini_freeway/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
To determine the weak solution to a non-convex scalar conservation law, in this problem, we will use Olenik entropy condition:
$$ \frac{f(\rho)-f(\rho_l)}{\rho-\rho_l}\geq \frac{f(\rho_l)-f(\rho_r)}{\rho_l-\rho_r}\geq \frac{f(\rho)-f(\rho_r)}{\rho-\rho_r} $$for all $\rho$ between $\rho_l$ and $\rho_r$.
The entropy-satisfying solution to this problem can be determined form the graph of $f_\epsilon(\rho)$. This method is called convex hull.
We use the 2nd order Lax-Wendroff method with Van Leer limiter.
We only focus on each car $x_n(t)$.
The density between two car is the reciprocal of the distance: $$\rho_{n+1}(t)=\frac{1}{x_{n}(t)-x_{n+1}(t)};$$
The velocity of each car only depends on the distance: $$v_{n}(t)=v(\rho_{n}(t));$$
The position is updated by: $$x_{n}(t+\Delta t)=x_{n}(t)+v_{n}(t)\Delta t.$$
Based on these assumptions, we can build up our own ODE Solver.
# ODE Solver Initialization
def initialization(N,rhol,rhor):
return linspace(0.,-(N-1.)/rhol,N)
# Main Function of ODE Solver
def ode_solver_freeway(N,tfinal,output_number,rhol,rhor):
from copy import copy
Dt = tfinal/output_number
dt = min(0.01,Dt)
x = zeros((output_number+1,N))
rho = copy(x)
x[0] = initialization(N,rhol,rhor)
rho[0][1:N+1] = -1./diff(x[0])
rho[0][0] = rhor
m = int(Dt/dt)
xs ,rhos = copy(x[0]), copy(rho[0])
for i in range(1,output_number+1):
for j in range(m):
v = flux_freeway_origin(rhos)/rhos
xs = xs + v*dt
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
v = flux_freeway_origin(rhos)/rhos
xs = xs + v*(Dt-dt*m)
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
x[i], rho[i] = copy(xs), copy(rhos)
return x, rho
N = 80
tfinal = 40.
output_number = 80
We construct the smallest convex hull of the set $\{(\rho,y):\rho_r<\rho<\rho_l,\ \text{and}\ y\leq f_\epsilon(\rho)\}$. When $\epsilon\rightarrow 0$, the rarefaction wave flattens out and reduces to the constant state $\rho_m$. The limiting solution is:
$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_m \ \ \ \text{if}\ \ st \leq x \leq t,\\ \rho_r \ \ \ \text{if}\ \ x > t, \end{array} \right.$$where $s=\frac{f_\epsilon(\rho_\star)-f_\epsilon(\rho_l)}{\rho_\star-\rho_l}=f'_\epsilon(\rho_{\star})$, and $\rho\in(\rho_m-\epsilon,\rho_m)$, since the shock and characteristic speeds are equal.
In the numerical simulation, the initial condition we are using is:
$$\rho_l=0.9, \rho_r=0.2.$$
%%system
make data ql=0.9 qr=0.2
make .plots
['rm -f .data', 'Wrote data to initial_data.txt', 'python setrun.py classic ', 'touch .data', 'gfortran -c qinit.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o qinit.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c setprob.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o setprob.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c zeroin.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o zeroin.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c rp1_freeway.f90 \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o rp1_freeway.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/setaux.f90 \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/setaux.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/bc1.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/bc1.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/b4step1.f90 \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/b4step1.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/driver.f90 \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/driver.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/claw1ez.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/claw1ez.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/claw1.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/claw1.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/copyq1.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/copyq1.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/inlinelimiter.f90 \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/inlinelimiter.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/opendatafile.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/opendatafile.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/out1.f \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/out1.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/src1.f90 \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/src1.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran -c /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/step1.f90 \t\t\t\t\t -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/step1.o', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', 'gfortran qinit.o setprob.o zeroin.o rp1_freeway.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/setaux.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/bc1.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/b4step1.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/driver.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/claw1ez.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/claw1.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/copyq1.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/inlinelimiter.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/opendatafile.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/out1.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/src1.o /Users/Peng/Documents/AMATH574/clawpack/classic/src/1d/step1.o -I/Users/Peng/Documents/AMATH574/am574-group07/presentation/ -o xclaw ', 'gfortran: warning: couldn\xe2\x80\x99t understand kern.osversion \xe2\x80\x9814.1.0', '/Library/Developer/CommandLineTools/usr/bin/make output', 'rm -f .output', 'python /Users/Peng/Documents/AMATH574/clawpack/clawutil/src/python/clawutil/runclaw.py xclaw _output \\', '\tTrue False . False False None', 'Reading data file: claw.data ', ' first 5 lines are comments and will be skipped', ' running...', ' ', 'Reading data file: setprob.data', ' first 5 lines are comments and will be skipped', 'CLAW1EZ: Frame 1 output files done at time t = 0.5000D-01', '', 'CLAW1EZ: Frame 2 output files done at time t = 0.1000D+00', '', 'CLAW1EZ: Frame 3 output files done at time t = 0.1500D+00', '', 'CLAW1EZ: Frame 4 output files done at time t = 0.2000D+00', '', 'CLAW1EZ: Frame 5 output files done at time t = 0.2500D+00', '', 'CLAW1EZ: Frame 6 output files done at time t = 0.3000D+00', '', 'CLAW1EZ: Frame 7 output files done at time t = 0.3500D+00', '', 'CLAW1EZ: Frame 8 output files done at time t = 0.4000D+00', '', 'CLAW1EZ: Frame 9 output files done at time t = 0.4500D+00', '', 'CLAW1EZ: Frame 10 output files done at time t = 0.5000D+00', '', 'CLAW1EZ: Frame 11 output files done at time t = 0.5500D+00', '', 'CLAW1EZ: Frame 12 output files done at time t = 0.6000D+00', '', 'CLAW1EZ: Frame 13 output files done at time t = 0.6500D+00', '', 'CLAW1EZ: Frame 14 output files done at time t = 0.7000D+00', '', 'CLAW1EZ: Frame 15 output files done at time t = 0.7500D+00', '', 'CLAW1EZ: Frame 16 output files done at time t = 0.8000D+00', '', 'CLAW1EZ: Frame 17 output files done at time t = 0.8500D+00', '', 'CLAW1EZ: Frame 18 output files done at time t = 0.9000D+00', '', 'CLAW1EZ: Frame 19 output files done at time t = 0.9500D+00', '', 'CLAW1EZ: Frame 20 output files done at time t = 0.1000D+01', '', 'CLAW1EZ: Frame 21 output files done at time t = 0.1050D+01', '', 'CLAW1EZ: Frame 22 output files done at time t = 0.1100D+01', '', 'CLAW1EZ: Frame 23 output files done at time t = 0.1150D+01', '', 'CLAW1EZ: Frame 24 output files done at time t = 0.1200D+01', '', 'CLAW1EZ: Frame 25 output files done at time t = 0.1250D+01', '', 'CLAW1EZ: Frame 26 output files done at time t = 0.1300D+01', '', 'CLAW1EZ: Frame 27 output files done at time t = 0.1350D+01', '', 'CLAW1EZ: Frame 28 output files done at time t = 0.1400D+01', '', 'CLAW1EZ: Frame 29 output files done at time t = 0.1450D+01', '', 'CLAW1EZ: Frame 30 output files done at time t = 0.1500D+01', '', 'CLAW1EZ: Frame 31 output files done at time t = 0.1550D+01', '', 'CLAW1EZ: Frame 32 output files done at time t = 0.1600D+01', '', 'CLAW1EZ: Frame 33 output files done at time t = 0.1650D+01', '', 'CLAW1EZ: Frame 34 output files done at time t = 0.1700D+01', '', 'CLAW1EZ: Frame 35 output files done at time t = 0.1750D+01', '', 'CLAW1EZ: Frame 36 output files done at time t = 0.1800D+01', '', 'CLAW1EZ: Frame 37 output files done at time t = 0.1850D+01', '', 'CLAW1EZ: Frame 38 output files done at time t = 0.1900D+01', '', 'CLAW1EZ: Frame 39 output files done at time t = 0.1950D+01', '', 'CLAW1EZ: Frame 40 output files done at time t = 0.2000D+01', '', 'Note: The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL', '==> runclaw: Will take data from /Users/Peng/Documents/AMATH574/am574-group07/presentation', '==> runclaw: Will write output to /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '==> runclaw: Removing all old fort files in /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '', '==> Running with command:', ' /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw', '', '==> runclaw: Finished executing', '', '==> runclaw: Done executing /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw via clawutil.runclaw.py', '==> runclaw: Output is in /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '/Library/Developer/CommandLineTools/usr/bin/make plots', 'rm -f .plots', 'python /Users/Peng/Documents/AMATH574/clawpack/visclaw/src/python/visclaw/plotclaw.py _output _plots setplot.py ', 'Importing setplot.setplot from /Users/Peng/Documents/AMATH574/am574-group07/presentation.', 'Executed setplot successfully', 'Will plot 41 frames numbered: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40]', 'Will make 1 figure(s) for each frame, numbered: [1]', '', '-----------------------------------', '', '', 'Creating html pages for figures...', '', "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ", ' already exists, files may be overwritten ', 'Now making png files for all figures...', ' Reading Frame 0 at t = 0 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 0 at time t = 0.0', ' Reading Frame 1 at t = 0.05 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 1 at time t = 0.05', ' Reading Frame 2 at t = 0.1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 2 at time t = 0.1', ' Reading Frame 3 at t = 0.15 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 3 at time t = 0.15', ' Reading Frame 4 at t = 0.2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 4 at time t = 0.2', ' Reading Frame 5 at t = 0.25 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 5 at time t = 0.25', ' Reading Frame 6 at t = 0.3 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 6 at time t = 0.3', ' Reading Frame 7 at t = 0.35 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 7 at time t = 0.35', ' Reading Frame 8 at t = 0.4 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 8 at time t = 0.4', ' Reading Frame 9 at t = 0.45 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 9 at time t = 0.45', ' Reading Frame 10 at t = 0.5 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 10 at time t = 0.5', ' Reading Frame 11 at t = 0.55 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 11 at time t = 0.55', ' Reading Frame 12 at t = 0.6 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 12 at time t = 0.6', ' Reading Frame 13 at t = 0.65 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 13 at time t = 0.65', ' Reading Frame 14 at t = 0.7 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 14 at time t = 0.7', ' Reading Frame 15 at t = 0.75 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 15 at time t = 0.75', ' Reading Frame 16 at t = 0.8 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 16 at time t = 0.8', ' Reading Frame 17 at t = 0.85 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 17 at time t = 0.85', ' Reading Frame 18 at t = 0.9 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 18 at time t = 0.9', ' Reading Frame 19 at t = 0.95 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 19 at time t = 0.95', ' Reading Frame 20 at t = 1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 20 at time t = 1.0', ' Reading Frame 21 at t = 1.05 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 21 at time t = 1.05', ' Reading Frame 22 at t = 1.1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 22 at time t = 1.1', ' Reading Frame 23 at t = 1.15 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 23 at time t = 1.15', ' Reading Frame 24 at t = 1.2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 24 at time t = 1.2', ' Reading Frame 25 at t = 1.25 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 25 at time t = 1.25', ' Reading Frame 26 at t = 1.3 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 26 at time t = 1.3', ' Reading Frame 27 at t = 1.35 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 27 at time t = 1.35', ' Reading Frame 28 at t = 1.4 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 28 at time t = 1.4', ' Reading Frame 29 at t = 1.45 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 29 at time t = 1.45', ' Reading Frame 30 at t = 1.5 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 30 at time t = 1.5', ' Reading Frame 31 at t = 1.55 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 31 at time t = 1.55', ' Reading Frame 32 at t = 1.6 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 32 at time t = 1.6', ' Reading Frame 33 at t = 1.65 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 33 at time t = 1.65', ' Reading Frame 34 at t = 1.7 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 34 at time t = 1.7', ' Reading Frame 35 at t = 1.75 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 35 at time t = 1.75', ' Reading Frame 36 at t = 1.8 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 36 at time t = 1.8', ' Reading Frame 37 at t = 1.85 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 37 at time t = 1.85', ' Reading Frame 38 at t = 1.9 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 38 at time t = 1.9', ' Reading Frame 39 at t = 1.95 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 39 at time t = 1.95', ' Reading Frame 40 at t = 2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 40 at time t = 2.0', '', '-----------------------------------', '', 'Creating latex file...', "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ", ' already exists, files may be overwritten ', '', 'Latex file created: ', ' /Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/plots.tex', '', 'Use pdflatex to create pdf file', 'Created JSAnimation for figure 1', '', '--------------------------------------------------------', '', 'Point your browser to:', ' file:///Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/_PlotIndex.html']
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
rhol = 0.9
rhor = 0.2
x, rho = ode_solver_freeway(N,tfinal,output_number,rhol,rhor)
%%system
mkdir _plots_ode_freeway1
['mkdir: _plots_ode_freeway1: File exists']
Dt = tfinal/output_number
for i in range(output_number+1):
clf()
figure(figsize=(14,5))
subplot(1,2,1)
for j in range(0,N,4):
plot(x.T[j][0:i+1],linspace(0.,Dt*i,i+1),'b');
plot(linspace(1./rhor,1./rhor+Dt*i,i+1),linspace(0.,Dt*i,i+1),'b')
plot(linspace(2./rhor,2./rhor+Dt*i,i+1),linspace(0.,Dt*i,i+1),'b')
plot(hstack((x[i][0:N:4],1./rhor+Dt*i,2./rhor+Dt*i)),Dt*i*ones(N/4+2),'ob')
plot([-N/rhol-10,tfinal+20],[Dt*i,Dt*i],'--r');
axis([-N/rhol-10,tfinal+20,0,tfinal+1])
title('Trajectory of Cars')
xlabel('$x$')
ylabel('$t$')
subplot(1,2,2)
plot(hstack((tfinal+10,x[i],-N/rhol-10)),hstack((rhor,rho[i],rhol)));
axis([-N/rhol-10,tfinal+10,-0.1,1.1])
title('Density')
xlabel('$x$')
ylabel('rho')
fname = '_plots_ode_freeway1/frame%04i'%i + 'fig1' + '.png'
savefig(fname)
<matplotlib.figure.Figure at 0x10ba35350>
<matplotlib.figure.Figure at 0x10ba35090>
<matplotlib.figure.Figure at 0x10bacc810>
<matplotlib.figure.Figure at 0x110714690>
<matplotlib.figure.Figure at 0x112cf29d0>
<matplotlib.figure.Figure at 0x112d90fd0>
<matplotlib.figure.Figure at 0x10b367d50>
<matplotlib.figure.Figure at 0x10d49c290>
<matplotlib.figure.Figure at 0x10bc0c190>
<matplotlib.figure.Figure at 0x10b46b3d0>
<matplotlib.figure.Figure at 0x10cf86c90>
<matplotlib.figure.Figure at 0x10badc850>
<matplotlib.figure.Figure at 0x10a3776d0>
<matplotlib.figure.Figure at 0x10b3cc710>
<matplotlib.figure.Figure at 0x10b6da110>
<matplotlib.figure.Figure at 0x10b3a3d90>
<matplotlib.figure.Figure at 0x106d79b90>
<matplotlib.figure.Figure at 0x10b8c2650>
<matplotlib.figure.Figure at 0x10b3d7e50>
<matplotlib.figure.Figure at 0x112cb0f50>
<matplotlib.figure.Figure at 0x10bb42510>
<matplotlib.figure.Figure at 0x10babfb10>
<matplotlib.figure.Figure at 0x10bc0c790>
<matplotlib.figure.Figure at 0x10b6b1290>
<matplotlib.figure.Figure at 0x10bb71d90>
<matplotlib.figure.Figure at 0x112c85510>
<matplotlib.figure.Figure at 0x1081e2790>
<matplotlib.figure.Figure at 0x10a3773d0>
<matplotlib.figure.Figure at 0x10bbd5150>
<matplotlib.figure.Figure at 0x10b404d90>
<matplotlib.figure.Figure at 0x10b8f5790>
<matplotlib.figure.Figure at 0x117f18f90>
<matplotlib.figure.Figure at 0x10b46b710>
<matplotlib.figure.Figure at 0x10b3d78d0>
<matplotlib.figure.Figure at 0x10cf792d0>
<matplotlib.figure.Figure at 0x10b6719d0>
<matplotlib.figure.Figure at 0x112d2c750>
<matplotlib.figure.Figure at 0x10bbd5250>
<matplotlib.figure.Figure at 0x106e4e410>
<matplotlib.figure.Figure at 0x10b3d7490>
<matplotlib.figure.Figure at 0x11419dc50>
<matplotlib.figure.Figure at 0x10bc4c990>
<matplotlib.figure.Figure at 0x10a351d50>
<matplotlib.figure.Figure at 0x106ce15d0>
<matplotlib.figure.Figure at 0x10b6a6d50>
<matplotlib.figure.Figure at 0x106c6a6d0>
<matplotlib.figure.Figure at 0x10b420950>
<matplotlib.figure.Figure at 0x10d49c650>
<matplotlib.figure.Figure at 0x112d06690>
<matplotlib.figure.Figure at 0x10b5abd10>
<matplotlib.figure.Figure at 0x10b0e5290>
<matplotlib.figure.Figure at 0x10b8b77d0>
<matplotlib.figure.Figure at 0x10bbac150>
<matplotlib.figure.Figure at 0x10822e350>
<matplotlib.figure.Figure at 0x10b6b1e90>
<matplotlib.figure.Figure at 0x112d26d50>
<matplotlib.figure.Figure at 0x10b61e310>
<matplotlib.figure.Figure at 0x112d93710>
<matplotlib.figure.Figure at 0x112cb8710>
<matplotlib.figure.Figure at 0x10b337fd0>
<matplotlib.figure.Figure at 0x10b307ad0>
<matplotlib.figure.Figure at 0x106d32690>
<matplotlib.figure.Figure at 0x10babf8d0>
<matplotlib.figure.Figure at 0x1121fa890>
<matplotlib.figure.Figure at 0x114101910>
<matplotlib.figure.Figure at 0x112d06990>
<matplotlib.figure.Figure at 0x10cf8b190>
<matplotlib.figure.Figure at 0x1082426d0>
<matplotlib.figure.Figure at 0x1164c0510>
<matplotlib.figure.Figure at 0x11075f210>
<matplotlib.figure.Figure at 0x10c1685d0>
<matplotlib.figure.Figure at 0x10b8c7390>
<matplotlib.figure.Figure at 0x1144c2c50>
<matplotlib.figure.Figure at 0x10f11cfd0>
<matplotlib.figure.Figure at 0x10e787210>
<matplotlib.figure.Figure at 0x10baa6190>
<matplotlib.figure.Figure at 0x10bc530d0>
<matplotlib.figure.Figure at 0x117f04ad0>
<matplotlib.figure.Figure at 0x114111350>
<matplotlib.figure.Figure at 0x107101f50>
<matplotlib.figure.Figure at 0x107765b90>
figno = 1
fname = '_plots_ode_freeway1/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(14, 5),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
We construct the smallest convex hull of the set $\{(\rho,y):\rho_l<\rho<\rho_r,\ \text{and}\ y\geq f_\epsilon(\rho)\}$. When $\epsilon\rightarrow 0$, the rarefaction wave flattens out and reduces to the constant state $\rho_m$. The limiting solution is:
$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_m \ \ \ \text{if}\ \ st \leq x \leq -\gamma t,\\ \rho_r \ \ \ \text{if}\ \ x > -\gamma t, \end{array} \right.$$where $s=\frac{f_\epsilon(\rho_\star)-f_\epsilon(\rho_l)}{\rho_\star-\rho_l}=f'_\epsilon(\rho_{\star})$, and $\rho\in(\rho_m-\epsilon,\rho_m)$.
In the numerical simulation, the initial condition we are using is:
$$\rho_l=0.4, \rho_r=0.9.$$
%%system
make data ql=0.4 qr=0.9
make .plots
['rm -f .data', 'Wrote data to initial_data.txt', 'python setrun.py classic ', 'touch .data', '/Library/Developer/CommandLineTools/usr/bin/make output', 'rm -f .output', 'python /Users/Peng/Documents/AMATH574/clawpack/clawutil/src/python/clawutil/runclaw.py xclaw _output \\', '\tTrue False . False False None', 'Reading data file: claw.data ', ' first 5 lines are comments and will be skipped', ' running...', ' ', 'Reading data file: setprob.data', ' first 5 lines are comments and will be skipped', 'CLAW1EZ: Frame 1 output files done at time t = 0.5000D-01', '', 'CLAW1EZ: Frame 2 output files done at time t = 0.1000D+00', '', 'CLAW1EZ: Frame 3 output files done at time t = 0.1500D+00', '', 'CLAW1EZ: Frame 4 output files done at time t = 0.2000D+00', '', 'CLAW1EZ: Frame 5 output files done at time t = 0.2500D+00', '', 'CLAW1EZ: Frame 6 output files done at time t = 0.3000D+00', '', 'CLAW1EZ: Frame 7 output files done at time t = 0.3500D+00', '', 'CLAW1EZ: Frame 8 output files done at time t = 0.4000D+00', '', 'CLAW1EZ: Frame 9 output files done at time t = 0.4500D+00', '', 'CLAW1EZ: Frame 10 output files done at time t = 0.5000D+00', '', 'CLAW1EZ: Frame 11 output files done at time t = 0.5500D+00', '', 'CLAW1EZ: Frame 12 output files done at time t = 0.6000D+00', '', 'CLAW1EZ: Frame 13 output files done at time t = 0.6500D+00', '', 'CLAW1EZ: Frame 14 output files done at time t = 0.7000D+00', '', 'CLAW1EZ: Frame 15 output files done at time t = 0.7500D+00', '', 'CLAW1EZ: Frame 16 output files done at time t = 0.8000D+00', '', 'CLAW1EZ: Frame 17 output files done at time t = 0.8500D+00', '', 'CLAW1EZ: Frame 18 output files done at time t = 0.9000D+00', '', 'CLAW1EZ: Frame 19 output files done at time t = 0.9500D+00', '', 'CLAW1EZ: Frame 20 output files done at time t = 0.1000D+01', '', 'CLAW1EZ: Frame 21 output files done at time t = 0.1050D+01', '', 'CLAW1EZ: Frame 22 output files done at time t = 0.1100D+01', '', 'CLAW1EZ: Frame 23 output files done at time t = 0.1150D+01', '', 'CLAW1EZ: Frame 24 output files done at time t = 0.1200D+01', '', 'CLAW1EZ: Frame 25 output files done at time t = 0.1250D+01', '', 'CLAW1EZ: Frame 26 output files done at time t = 0.1300D+01', '', 'CLAW1EZ: Frame 27 output files done at time t = 0.1350D+01', '', 'CLAW1EZ: Frame 28 output files done at time t = 0.1400D+01', '', 'CLAW1EZ: Frame 29 output files done at time t = 0.1450D+01', '', 'CLAW1EZ: Frame 30 output files done at time t = 0.1500D+01', '', 'CLAW1EZ: Frame 31 output files done at time t = 0.1550D+01', '', 'CLAW1EZ: Frame 32 output files done at time t = 0.1600D+01', '', 'CLAW1EZ: Frame 33 output files done at time t = 0.1650D+01', '', 'CLAW1EZ: Frame 34 output files done at time t = 0.1700D+01', '', 'CLAW1EZ: Frame 35 output files done at time t = 0.1750D+01', '', 'CLAW1EZ: Frame 36 output files done at time t = 0.1800D+01', '', 'CLAW1EZ: Frame 37 output files done at time t = 0.1850D+01', '', 'CLAW1EZ: Frame 38 output files done at time t = 0.1900D+01', '', 'CLAW1EZ: Frame 39 output files done at time t = 0.1950D+01', '', 'CLAW1EZ: Frame 40 output files done at time t = 0.2000D+01', '', 'Note: The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL', '==> runclaw: Will take data from /Users/Peng/Documents/AMATH574/am574-group07/presentation', '==> runclaw: Will write output to /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '==> runclaw: Removing all old fort files in /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '', '==> Running with command:', ' /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw', '', '==> runclaw: Finished executing', '', '==> runclaw: Done executing /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw via clawutil.runclaw.py', '==> runclaw: Output is in /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '/Library/Developer/CommandLineTools/usr/bin/make plots', 'rm -f .plots', 'python /Users/Peng/Documents/AMATH574/clawpack/visclaw/src/python/visclaw/plotclaw.py _output _plots setplot.py ', 'Importing setplot.setplot from /Users/Peng/Documents/AMATH574/am574-group07/presentation.', 'Executed setplot successfully', 'Will plot 41 frames numbered: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40]', 'Will make 1 figure(s) for each frame, numbered: [1]', '', '-----------------------------------', '', '', 'Creating html pages for figures...', '', "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ", ' already exists, files may be overwritten ', 'Now making png files for all figures...', ' Reading Frame 0 at t = 0 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 0 at time t = 0.0', ' Reading Frame 1 at t = 0.05 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 1 at time t = 0.05', ' Reading Frame 2 at t = 0.1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 2 at time t = 0.1', ' Reading Frame 3 at t = 0.15 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 3 at time t = 0.15', ' Reading Frame 4 at t = 0.2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 4 at time t = 0.2', ' Reading Frame 5 at t = 0.25 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 5 at time t = 0.25', ' Reading Frame 6 at t = 0.3 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 6 at time t = 0.3', ' Reading Frame 7 at t = 0.35 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 7 at time t = 0.35', ' Reading Frame 8 at t = 0.4 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 8 at time t = 0.4', ' Reading Frame 9 at t = 0.45 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 9 at time t = 0.45', ' Reading Frame 10 at t = 0.5 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 10 at time t = 0.5', ' Reading Frame 11 at t = 0.55 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 11 at time t = 0.55', ' Reading Frame 12 at t = 0.6 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 12 at time t = 0.6', ' Reading Frame 13 at t = 0.65 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 13 at time t = 0.65', ' Reading Frame 14 at t = 0.7 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 14 at time t = 0.7', ' Reading Frame 15 at t = 0.75 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 15 at time t = 0.75', ' Reading Frame 16 at t = 0.8 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 16 at time t = 0.8', ' Reading Frame 17 at t = 0.85 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 17 at time t = 0.85', ' Reading Frame 18 at t = 0.9 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 18 at time t = 0.9', ' Reading Frame 19 at t = 0.95 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 19 at time t = 0.95', ' Reading Frame 20 at t = 1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 20 at time t = 1.0', ' Reading Frame 21 at t = 1.05 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 21 at time t = 1.05', ' Reading Frame 22 at t = 1.1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 22 at time t = 1.1', ' Reading Frame 23 at t = 1.15 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 23 at time t = 1.15', ' Reading Frame 24 at t = 1.2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 24 at time t = 1.2', ' Reading Frame 25 at t = 1.25 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 25 at time t = 1.25', ' Reading Frame 26 at t = 1.3 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 26 at time t = 1.3', ' Reading Frame 27 at t = 1.35 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 27 at time t = 1.35', ' Reading Frame 28 at t = 1.4 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 28 at time t = 1.4', ' Reading Frame 29 at t = 1.45 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 29 at time t = 1.45', ' Reading Frame 30 at t = 1.5 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 30 at time t = 1.5', ' Reading Frame 31 at t = 1.55 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 31 at time t = 1.55', ' Reading Frame 32 at t = 1.6 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 32 at time t = 1.6', ' Reading Frame 33 at t = 1.65 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 33 at time t = 1.65', ' Reading Frame 34 at t = 1.7 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 34 at time t = 1.7', ' Reading Frame 35 at t = 1.75 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 35 at time t = 1.75', ' Reading Frame 36 at t = 1.8 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 36 at time t = 1.8', ' Reading Frame 37 at t = 1.85 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 37 at time t = 1.85', ' Reading Frame 38 at t = 1.9 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 38 at time t = 1.9', ' Reading Frame 39 at t = 1.95 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 39 at time t = 1.95', ' Reading Frame 40 at t = 2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 40 at time t = 2.0', '', '-----------------------------------', '', 'Creating latex file...', "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ", ' already exists, files may be overwritten ', '', 'Latex file created: ', ' /Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/plots.tex', '', 'Use pdflatex to create pdf file', 'Created JSAnimation for figure 1', '', '--------------------------------------------------------', '', 'Point your browser to:', ' file:///Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/_PlotIndex.html']
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
rhol = 0.4
rhor = 0.9
x, rho = ode_solver_freeway(N,tfinal,output_number,rhol,rhor)
%%system
mkdir _plots_ode_freeway2
['mkdir: _plots_ode_freeway2: File exists']
Dt = tfinal/output_number
for i in range(output_number+1):
clf()
figure(figsize=(14,5))
subplot(1,2,1)
for j in range(0,N,4):
plot(x.T[j][0:i+1],linspace(0.,Dt*i,i+1),'b');
v = gamma*(1.-rhor)/rhor
plot(linspace(1./rhor,1./rhor+Dt*i*v,i+1),linspace(0.,Dt*i,i+1),'b')
plot(linspace(2./rhor,2./rhor+Dt*i*v,i+1),linspace(0.,Dt*i,i+1),'b')
plot(hstack((x[i][0:N:4],1./rhor+Dt*i*v,2./rhor+Dt*i*v)),Dt*i*ones(N/4+2),'ob')
plot([-N/rhol-10,tfinal+20],[Dt*i,Dt*i],'--r');
axis([-N/rhol-10,tfinal+20,0,tfinal+1])
title('Trajectory of Cars')
xlabel('$x$')
ylabel('$t$')
subplot(1,2,2)
plot(hstack((tfinal+10,x[i],-N/rhol-10)),hstack((rhor,rho[i],rhol)));
axis([-N/rhol-10,tfinal+10,-0.1,1.1])
title('Density')
xlabel('$x$')
ylabel('rho')
fname = '_plots_ode_freeway2/frame%04i'%i + 'fig1' + '.png'
savefig(fname)
<matplotlib.figure.Figure at 0x111276850>
<matplotlib.figure.Figure at 0x111276c90>
<matplotlib.figure.Figure at 0x1074b2150>
<matplotlib.figure.Figure at 0x10bc53c10>
<matplotlib.figure.Figure at 0x10b46bb10>
<matplotlib.figure.Figure at 0x10bb498d0>
<matplotlib.figure.Figure at 0x10b8c7e50>
<matplotlib.figure.Figure at 0x10b469710>
<matplotlib.figure.Figure at 0x106fa8790>
<matplotlib.figure.Figure at 0x106f92250>
<matplotlib.figure.Figure at 0x106d87210>
<matplotlib.figure.Figure at 0x1112489d0>
<matplotlib.figure.Figure at 0x10a0e8c50>
<matplotlib.figure.Figure at 0x112d97d90>
<matplotlib.figure.Figure at 0x107626fd0>
<matplotlib.figure.Figure at 0x10777e550>
<matplotlib.figure.Figure at 0x10f11cf10>
<matplotlib.figure.Figure at 0x1164c9250>
<matplotlib.figure.Figure at 0x112d26490>
<matplotlib.figure.Figure at 0x11418f590>
<matplotlib.figure.Figure at 0x10ba3d190>
<matplotlib.figure.Figure at 0x10bb91950>
<matplotlib.figure.Figure at 0x10b3437d0>
<matplotlib.figure.Figure at 0x112c85dd0>
<matplotlib.figure.Figure at 0x1163ae110>
<matplotlib.figure.Figure at 0x10d4e4fd0>
<matplotlib.figure.Figure at 0x112cbaed0>
<matplotlib.figure.Figure at 0x10b46ef90>
<matplotlib.figure.Figure at 0x117f06e50>
<matplotlib.figure.Figure at 0x1081e2790>
<matplotlib.figure.Figure at 0x117ee39d0>
<matplotlib.figure.Figure at 0x110764f90>
<matplotlib.figure.Figure at 0x10b888410>
<matplotlib.figure.Figure at 0x11127c8d0>
<matplotlib.figure.Figure at 0x10b3bd6d0>
<matplotlib.figure.Figure at 0x112d82050>
<matplotlib.figure.Figure at 0x10b333890>
<matplotlib.figure.Figure at 0x1062daa50>
<matplotlib.figure.Figure at 0x114497990>
<matplotlib.figure.Figure at 0x111255b50>
<matplotlib.figure.Figure at 0x10b19bd10>
<matplotlib.figure.Figure at 0x10824b4d0>
<matplotlib.figure.Figure at 0x106ef9f50>
<matplotlib.figure.Figure at 0x10b3d7750>
<matplotlib.figure.Figure at 0x10b621890>
<matplotlib.figure.Figure at 0x10b8b76d0>
<matplotlib.figure.Figure at 0x10b6c31d0>
<matplotlib.figure.Figure at 0x108f026d0>
<matplotlib.figure.Figure at 0x10cfa6ed0>
<matplotlib.figure.Figure at 0x108240350>
<matplotlib.figure.Figure at 0x107332f90>
<matplotlib.figure.Figure at 0x1163ebd90>
<matplotlib.figure.Figure at 0x10bb01590>
<matplotlib.figure.Figure at 0x106f9d550>
<matplotlib.figure.Figure at 0x1141014d0>
<matplotlib.figure.Figure at 0x11649b590>
<matplotlib.figure.Figure at 0x11639f450>
<matplotlib.figure.Figure at 0x112d932d0>
<matplotlib.figure.Figure at 0x106f48450>
<matplotlib.figure.Figure at 0x107332410>
<matplotlib.figure.Figure at 0x10b110a90>
<matplotlib.figure.Figure at 0x111245f10>
<matplotlib.figure.Figure at 0x1074d1ed0>
<matplotlib.figure.Figure at 0x10b586690>
<matplotlib.figure.Figure at 0x10ba71090>
<matplotlib.figure.Figure at 0x106d79c50>
<matplotlib.figure.Figure at 0x10b62da10>
<matplotlib.figure.Figure at 0x10b848a90>
<matplotlib.figure.Figure at 0x10b390ad0>
<matplotlib.figure.Figure at 0x10bc40f50>
<matplotlib.figure.Figure at 0x10b586150>
<matplotlib.figure.Figure at 0x107410a50>
<matplotlib.figure.Figure at 0x12165edd0>
<matplotlib.figure.Figure at 0x111260990>
<matplotlib.figure.Figure at 0x10f155650>
<matplotlib.figure.Figure at 0x110756450>
<matplotlib.figure.Figure at 0x10e79edd0>
<matplotlib.figure.Figure at 0x10b81b950>
<matplotlib.figure.Figure at 0x108173e90>
<matplotlib.figure.Figure at 0x10b6cd090>
<matplotlib.figure.Figure at 0x110769a10>
figno = 1
fname = '_plots_ode_freeway2/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(14, 5),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
We construct the same convex hull as Case 2. There is only a shock in this solution.
$$\rho(x,t)=\left\{ \begin{array}{l l} \rho_l \ \ \ \text{if}\ \ x < st,\\ \rho_r \ \ \ \text{if}\ \ x > st, \end{array} \right.$$where $s=\frac{f_\epsilon(\rho_r)-f_\epsilon(\rho_l)}{\rho_r-\rho_l}.$
In the numerical simulation, the initial condition we are using is:
$$\rho_l=0.2, \rho_r=0.9.$$
%%system
make data ql=0.2 qr=0.9
make .plots
['rm -f .data', 'Wrote data to initial_data.txt', 'python setrun.py classic ', 'touch .data', '/Library/Developer/CommandLineTools/usr/bin/make output', 'rm -f .output', 'python /Users/Peng/Documents/AMATH574/clawpack/clawutil/src/python/clawutil/runclaw.py xclaw _output \\', '\tTrue False . False False None', 'Reading data file: claw.data ', ' first 5 lines are comments and will be skipped', ' running...', ' ', 'Reading data file: setprob.data', ' first 5 lines are comments and will be skipped', 'CLAW1EZ: Frame 1 output files done at time t = 0.5000D-01', '', 'CLAW1EZ: Frame 2 output files done at time t = 0.1000D+00', '', 'CLAW1EZ: Frame 3 output files done at time t = 0.1500D+00', '', 'CLAW1EZ: Frame 4 output files done at time t = 0.2000D+00', '', 'CLAW1EZ: Frame 5 output files done at time t = 0.2500D+00', '', 'CLAW1EZ: Frame 6 output files done at time t = 0.3000D+00', '', 'CLAW1EZ: Frame 7 output files done at time t = 0.3500D+00', '', 'CLAW1EZ: Frame 8 output files done at time t = 0.4000D+00', '', 'CLAW1EZ: Frame 9 output files done at time t = 0.4500D+00', '', 'CLAW1EZ: Frame 10 output files done at time t = 0.5000D+00', '', 'CLAW1EZ: Frame 11 output files done at time t = 0.5500D+00', '', 'CLAW1EZ: Frame 12 output files done at time t = 0.6000D+00', '', 'CLAW1EZ: Frame 13 output files done at time t = 0.6500D+00', '', 'CLAW1EZ: Frame 14 output files done at time t = 0.7000D+00', '', 'CLAW1EZ: Frame 15 output files done at time t = 0.7500D+00', '', 'CLAW1EZ: Frame 16 output files done at time t = 0.8000D+00', '', 'CLAW1EZ: Frame 17 output files done at time t = 0.8500D+00', '', 'CLAW1EZ: Frame 18 output files done at time t = 0.9000D+00', '', 'CLAW1EZ: Frame 19 output files done at time t = 0.9500D+00', '', 'CLAW1EZ: Frame 20 output files done at time t = 0.1000D+01', '', 'CLAW1EZ: Frame 21 output files done at time t = 0.1050D+01', '', 'CLAW1EZ: Frame 22 output files done at time t = 0.1100D+01', '', 'CLAW1EZ: Frame 23 output files done at time t = 0.1150D+01', '', 'CLAW1EZ: Frame 24 output files done at time t = 0.1200D+01', '', 'CLAW1EZ: Frame 25 output files done at time t = 0.1250D+01', '', 'CLAW1EZ: Frame 26 output files done at time t = 0.1300D+01', '', 'CLAW1EZ: Frame 27 output files done at time t = 0.1350D+01', '', 'CLAW1EZ: Frame 28 output files done at time t = 0.1400D+01', '', 'CLAW1EZ: Frame 29 output files done at time t = 0.1450D+01', '', 'CLAW1EZ: Frame 30 output files done at time t = 0.1500D+01', '', 'CLAW1EZ: Frame 31 output files done at time t = 0.1550D+01', '', 'CLAW1EZ: Frame 32 output files done at time t = 0.1600D+01', '', 'CLAW1EZ: Frame 33 output files done at time t = 0.1650D+01', '', 'CLAW1EZ: Frame 34 output files done at time t = 0.1700D+01', '', 'CLAW1EZ: Frame 35 output files done at time t = 0.1750D+01', '', 'CLAW1EZ: Frame 36 output files done at time t = 0.1800D+01', '', 'CLAW1EZ: Frame 37 output files done at time t = 0.1850D+01', '', 'CLAW1EZ: Frame 38 output files done at time t = 0.1900D+01', '', 'CLAW1EZ: Frame 39 output files done at time t = 0.1950D+01', '', 'CLAW1EZ: Frame 40 output files done at time t = 0.2000D+01', '', 'Note: The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL', '==> runclaw: Will take data from /Users/Peng/Documents/AMATH574/am574-group07/presentation', '==> runclaw: Will write output to /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '==> runclaw: Removing all old fort files in /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '', '==> Running with command:', ' /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw', '', '==> runclaw: Finished executing', '', '==> runclaw: Done executing /Users/Peng/Documents/AMATH574/am574-group07/presentation/xclaw via clawutil.runclaw.py', '==> runclaw: Output is in /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', '/Library/Developer/CommandLineTools/usr/bin/make plots', 'rm -f .plots', 'python /Users/Peng/Documents/AMATH574/clawpack/visclaw/src/python/visclaw/plotclaw.py _output _plots setplot.py ', 'Importing setplot.setplot from /Users/Peng/Documents/AMATH574/am574-group07/presentation.', 'Executed setplot successfully', 'Will plot 41 frames numbered: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40]', 'Will make 1 figure(s) for each frame, numbered: [1]', '', '-----------------------------------', '', '', 'Creating html pages for figures...', '', "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ", ' already exists, files may be overwritten ', 'Now making png files for all figures...', ' Reading Frame 0 at t = 0 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 0 at time t = 0.0', ' Reading Frame 1 at t = 0.05 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 1 at time t = 0.05', ' Reading Frame 2 at t = 0.1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 2 at time t = 0.1', ' Reading Frame 3 at t = 0.15 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 3 at time t = 0.15', ' Reading Frame 4 at t = 0.2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 4 at time t = 0.2', ' Reading Frame 5 at t = 0.25 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 5 at time t = 0.25', ' Reading Frame 6 at t = 0.3 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 6 at time t = 0.3', ' Reading Frame 7 at t = 0.35 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 7 at time t = 0.35', ' Reading Frame 8 at t = 0.4 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 8 at time t = 0.4', ' Reading Frame 9 at t = 0.45 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 9 at time t = 0.45', ' Reading Frame 10 at t = 0.5 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 10 at time t = 0.5', ' Reading Frame 11 at t = 0.55 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 11 at time t = 0.55', ' Reading Frame 12 at t = 0.6 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 12 at time t = 0.6', ' Reading Frame 13 at t = 0.65 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 13 at time t = 0.65', ' Reading Frame 14 at t = 0.7 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 14 at time t = 0.7', ' Reading Frame 15 at t = 0.75 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 15 at time t = 0.75', ' Reading Frame 16 at t = 0.8 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 16 at time t = 0.8', ' Reading Frame 17 at t = 0.85 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 17 at time t = 0.85', ' Reading Frame 18 at t = 0.9 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 18 at time t = 0.9', ' Reading Frame 19 at t = 0.95 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 19 at time t = 0.95', ' Reading Frame 20 at t = 1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 20 at time t = 1.0', ' Reading Frame 21 at t = 1.05 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 21 at time t = 1.05', ' Reading Frame 22 at t = 1.1 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 22 at time t = 1.1', ' Reading Frame 23 at t = 1.15 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 23 at time t = 1.15', ' Reading Frame 24 at t = 1.2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 24 at time t = 1.2', ' Reading Frame 25 at t = 1.25 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 25 at time t = 1.25', ' Reading Frame 26 at t = 1.3 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 26 at time t = 1.3', ' Reading Frame 27 at t = 1.35 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 27 at time t = 1.35', ' Reading Frame 28 at t = 1.4 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 28 at time t = 1.4', ' Reading Frame 29 at t = 1.45 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 29 at time t = 1.45', ' Reading Frame 30 at t = 1.5 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 30 at time t = 1.5', ' Reading Frame 31 at t = 1.55 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 31 at time t = 1.55', ' Reading Frame 32 at t = 1.6 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 32 at time t = 1.6', ' Reading Frame 33 at t = 1.65 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 33 at time t = 1.65', ' Reading Frame 34 at t = 1.7 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 34 at time t = 1.7', ' Reading Frame 35 at t = 1.75 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 35 at time t = 1.75', ' Reading Frame 36 at t = 1.8 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 36 at time t = 1.8', ' Reading Frame 37 at t = 1.85 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 37 at time t = 1.85', ' Reading Frame 38 at t = 1.9 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 38 at time t = 1.9', ' Reading Frame 39 at t = 1.95 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 39 at time t = 1.95', ' Reading Frame 40 at t = 2 from outdir = /Users/Peng/Documents/AMATH574/am574-group07/presentation/_output', 'Frame 40 at time t = 2.0', '', '-----------------------------------', '', 'Creating latex file...', "Directory '/Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots' ", ' already exists, files may be overwritten ', '', 'Latex file created: ', ' /Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/plots.tex', '', 'Use pdflatex to create pdf file', 'Created JSAnimation for figure 1', '', '--------------------------------------------------------', '', 'Point your browser to:', ' file:///Users/Peng/Documents/AMATH574/am574-group07/presentation/_plots/_PlotIndex.html']
figno = 1
fname = '_plots/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(10, 6),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
rhol = 0.2
rhor = 0.9
x, rho = ode_solver_freeway(N,tfinal,output_number,rhol,rhor)
%%system
mkdir _plots_ode_freeway3
['mkdir: _plots_ode_freeway3: File exists']
Dt = tfinal/output_number
for i in range(output_number+1):
clf()
figure(figsize=(14,5))
v = gamma*(1.-rhor)/rhor
subplot(1,2,1)
for j in range(0,N,4):
plot(x.T[j][0:i+1],linspace(0.,Dt*i,i+1),'b');
plot(linspace(1./rhor,1./rhor+Dt*i*v,i+1),linspace(0.,Dt*i,i+1),'b')
plot(linspace(2./rhor,2./rhor+Dt*i*v,i+1),linspace(0.,Dt*i,i+1),'b')
plot(hstack((x[i][0:N:4],1./rhor+Dt*i*v,2./rhor+Dt*i*v)),Dt*i*ones(N/4+2),'ob')
plot([-N/rhol-10,tfinal+20],[Dt*i,Dt*i],'--r');
axis([-N/rhol-10,tfinal+20,0,tfinal+1])
title('Trajectory of Cars')
xlabel('$x$')
ylabel('$t$')
subplot(1,2,2)
plot(hstack((tfinal+10,x[i],-N/rhol-10)),hstack((rhor,rho[i],rhol)));
axis([-N/rhol-10,tfinal+10,-0.1,1.1])
title('Density')
xlabel('$x$')
ylabel('rho')
fname = '_plots_ode_freeway3/frame%04i'%i + 'fig1' + '.png'
savefig(fname)
<matplotlib.figure.Figure at 0x10cf86350>
<matplotlib.figure.Figure at 0x10cf86890>
<matplotlib.figure.Figure at 0x10b68d910>
<matplotlib.figure.Figure at 0x11123a890>
<matplotlib.figure.Figure at 0x1112ae7d0>
<matplotlib.figure.Figure at 0x10b691590>
<matplotlib.figure.Figure at 0x10ba359d0>
<matplotlib.figure.Figure at 0x10bbc8f10>
<matplotlib.figure.Figure at 0x117f27610>
<matplotlib.figure.Figure at 0x106f48390>
<matplotlib.figure.Figure at 0x112cbaf10>
<matplotlib.figure.Figure at 0x10d4d3790>
<matplotlib.figure.Figure at 0x10bc0cfd0>
<matplotlib.figure.Figure at 0x10b870910>
<matplotlib.figure.Figure at 0x10b6be650>
<matplotlib.figure.Figure at 0x106f48550>
<matplotlib.figure.Figure at 0x10babf150>
<matplotlib.figure.Figure at 0x10bba4190>
<matplotlib.figure.Figure at 0x1164c9d90>
<matplotlib.figure.Figure at 0x1164aa490>
<matplotlib.figure.Figure at 0x10ba82990>
<matplotlib.figure.Figure at 0x112d88b90>
<matplotlib.figure.Figure at 0x10cf94ed0>
<matplotlib.figure.Figure at 0x107775e90>
<matplotlib.figure.Figure at 0x111264150>
<matplotlib.figure.Figure at 0x11128e850>
<matplotlib.figure.Figure at 0x10b19eed0>
<matplotlib.figure.Figure at 0x10a039710>
<matplotlib.figure.Figure at 0x11418fe50>
<matplotlib.figure.Figure at 0x1144cfa90>
<matplotlib.figure.Figure at 0x10815bc10>
<matplotlib.figure.Figure at 0x10b8c2310>
<matplotlib.figure.Figure at 0x10bb49e10>
<matplotlib.figure.Figure at 0x108eed150>
<matplotlib.figure.Figure at 0x117efc150>
<matplotlib.figure.Figure at 0x1141a10d0>
<matplotlib.figure.Figure at 0x11649b690>
<matplotlib.figure.Figure at 0x10bc4ca50>
<matplotlib.figure.Figure at 0x106eba450>
<matplotlib.figure.Figure at 0x10cf8de10>
<matplotlib.figure.Figure at 0x106e4e410>
<matplotlib.figure.Figure at 0x10ba2e350>
<matplotlib.figure.Figure at 0x117f04ed0>
<matplotlib.figure.Figure at 0x10cf98c90>
<matplotlib.figure.Figure at 0x10b414fd0>
<matplotlib.figure.Figure at 0x1163c7790>
<matplotlib.figure.Figure at 0x108f6ea50>
<matplotlib.figure.Figure at 0x12164e110>
<matplotlib.figure.Figure at 0x11125d350>
<matplotlib.figure.Figure at 0x10b8f5990>
<matplotlib.figure.Figure at 0x10ba82f90>
<matplotlib.figure.Figure at 0x10b322c50>
<matplotlib.figure.Figure at 0x10e77b9d0>
<matplotlib.figure.Figure at 0x10b3cc2d0>
<matplotlib.figure.Figure at 0x10b8b74d0>
<matplotlib.figure.Figure at 0x1071213d0>
<matplotlib.figure.Figure at 0x1216423d0>
<matplotlib.figure.Figure at 0x10b888810>
<matplotlib.figure.Figure at 0x10b687890>
<matplotlib.figure.Figure at 0x10b41b290>
<matplotlib.figure.Figure at 0x10b691750>
<matplotlib.figure.Figure at 0x10bb54b90>
<matplotlib.figure.Figure at 0x107765f10>
<matplotlib.figure.Figure at 0x10c1483d0>
<matplotlib.figure.Figure at 0x10b414150>
<matplotlib.figure.Figure at 0x1144cfb90>
<matplotlib.figure.Figure at 0x108f199d0>
<matplotlib.figure.Figure at 0x10ba82bd0>
<matplotlib.figure.Figure at 0x11182ee10>
<matplotlib.figure.Figure at 0x10f163650>
<matplotlib.figure.Figure at 0x10b6cdd50>
<matplotlib.figure.Figure at 0x111267890>
<matplotlib.figure.Figure at 0x1164c41d0>
<matplotlib.figure.Figure at 0x107765210>
<matplotlib.figure.Figure at 0x10b388550>
<matplotlib.figure.Figure at 0x1163ae910>
<matplotlib.figure.Figure at 0x1117ce810>
<matplotlib.figure.Figure at 0x1077ae210>
<matplotlib.figure.Figure at 0x106eba110>
<matplotlib.figure.Figure at 0x10b8e4dd0>
<matplotlib.figure.Figure at 0x10bbc88d0>
figno = 1
fname = '_plots_ode_freeway3/*fig' + str(figno) + '.png'
filenames=sorted(glob.glob(fname))
fig = plt.figure(figsize=(14, 5),dpi=80)
ax = fig.add_axes([0, 0, 1, 1])
ax.axis('off')
im = plt.imshow(image.imread(filenames[0]))
animation.FuncAnimation(fig, animate, init_func=init,
frames=len(filenames), interval=200, blit=True)
When driving at night on an unfamiliar mountain road, facing with the empty road will cause the speed limitation by the uncertainty of the road condition. However if there are other cars in front of the road, it will be much easier to drive faster, since their tail lights indicate how the road twists and turns.
The velocity function for night-time driving is
$$ v(\rho)=\left\{ \begin{array}{l l} v_0 \ \ \ \text{if}\ \ \rho<\rho_a,\\ c\rho \ \ \ \text{if}\ \ \rho_a\leq \rho\leq \rho_b,\\ v_1(1-\rho) \ \ \ \text{if}\ \ \rho>\rho_b, \end{array} \right. $$Hence,
$$ f(\rho)=\left\{ \begin{array}{l l} v_0\rho \ \ \ \text{if}\ \ \rho<\rho_a,\\ c\rho^2 \ \ \ \text{if}\ \ \rho_a\leq \rho\leq \rho_b,\\ v_1\rho(1-\rho) \ \ \ \text{if}\ \ \rho>\rho_b, \end{array} \right. $$def vel_night(rho):
rhoa, rhob = 0.1, 0.3
v0, vmax = 1., 3.
v1 = vmax/(1.-rhob)
c = (vmax-v0)/(rhob-rhoa)
n = rho.size
v = zeros(n)
for i in range(n):
if rho[i] >= 0. and rho[i] < rhoa:
v[i] = v0
elif rho[i] >=rhoa and rho[i] <= rhob:
v[i] = c*rho[i]
elif rho[i] > rhob and rho[i] <= 1.:
v[i] = v1*(1.-rho[i])
else:
print "Invalid density value!"
return v
In the car-following model, the parameters are: $$ \rho_a=0.1,\rho_b=0.3, c=10, v_0=1. $$
rhopic = linspace(0.,1.,201)
figure(figsize=(14,6))
subplot(1,2,1)
plot(rhopic,vel_night(rhopic),'b')
xlabel('rho')
ylabel('$v$')
title('Velocity')
subplot(1,2,2)
plot(rhopic,rhopic*vel_night(rhopic),'b')
xlabel('rho')
ylabel('$f$')
title('Flux')
<matplotlib.text.Text at 0x106c861d0>
We take uniformly spaced cars as initial density on an otherwise empty road.
In LeVeque's paper:
"The lead driver sees $\rho=0$ and drives at speed $v_0$, and then the second driver initially sees $\rho=\rho_0$ and can drive faster. But as the second car accelerates, the third driver observes a drop in density and hence this driver starts to slow down. This causes the fourth driver to go faster."
It has been proven that the uniformly spaced traffic with distance $\frac{1}{\rho_0}$ with speed $v(\rho_0)$ is unstable under pertubations.
def ode_solver_night(N,tfinal,output_number,rhol,rhor):
from copy import copy
Dt = tfinal/output_number
dt = min(0.01,Dt)
x = zeros((output_number+1,N))
rho = copy(x)
x[0] = initialization(N,rhol,rhor)
rho[0][1:N+1] = -1./diff(x[0])
rho[0][0] = rhor
m = int(Dt/dt)
xs ,rhos = copy(x[0]), copy(rho[0])
for i in range(1,output_number+1):
for j in range(m):
v = vel_night(rhos)
xs = xs + v*dt
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
v = vel_night(rhos)
xs = xs + v*(Dt-dt*m)
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
x[i], rho[i] = copy(xs), copy(rhos)
return x, rho
N = 48
tfinal = 80
output_number = 80
figure(figsize=(15,5))
rhol = 1./12
rhor = 0.
subplot(1,3,1)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/12$')
rhol = 1./5
rhor = 0.
subplot(1,3,2)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/5$')
rhol = 1./3
rhor = 0.
subplot(1,3,3)
x, rho = ode_solver_night(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/3$')
<matplotlib.text.Text at 0x1074100d0>
With random perturbations to the speeds, 2-car platoons can no longer be observed. Instead, there are random n-car platoons, which is more reasonable to model the night-time traffic.
def ode_solver_night_rand(N,tfinal,output_number,rhol,rhor):
from copy import copy
from numpy.random import rand
A = 0.1
Dt = tfinal/output_number
dt = 0.01
x = zeros((output_number+1,N))
rho = copy(x)
x[0] = initialization(N,rhol,rhor)
rho[0][1:N+1] = -1./diff(x[0])
rho[0][0] = rhor
m = int(Dt/dt)
xs ,rhos = copy(x[0]), copy(rho[0])
for i in range(1,output_number+1):
for j in range(m):
v = vel_night(rhos) + A*(rand(N)-0.5*ones(N))
xs = xs + v*dt
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
v = vel_night(rhos) + A*(rand(N)-0.5*ones(N))
xs = xs + v*(Dt-dt*m)
rhos[1:N+1] = -1./diff(xs)
rhos[0] = rhor
x[i], rho[i] = copy(xs), copy(rhos)
return x, rho
N = 48
tfinal = 80
output_number = 80
figure(figsize=(15,5))
rhol = 1./12
rhor = 0.
subplot(1,3,1)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/12$')
rhol = 1./5
rhor = 0.
subplot(1,3,2)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/5$')
rhol = 1./3
rhor = 0.
subplot(1,3,3)
x, rho = ode_solver_night_rand(N,tfinal,output_number,rhol,rhor)
for i in range(N):
plot(x.T[i],linspace(0.,tfinal,output_number+1),'b')
axis([-200.,100.,0.,80.])
title('rho$=1/3$')
<matplotlib.text.Text at 0x10c18c3d0>